classdef instance_opt2   
    
    properties
        past_window_size=60;
        future_window_size=60;
        samplingTime = 1;
        system_=[];
        control_interval = 15; 
    end
    
    methods   
        function obj=instance_opt2

        end
            
         function [m_x,d_x]=linearize(obj,x)
            m_x = mean(x);  
            % d_x = x-m_x;  
           d_x = x-ones(size(x,1),1)*m_x;
        end
        
        function d_x=linearize_around(obj,x,m_x)
            d_x = x-ones( size(x,1),1) * m_x ;
        end
        
        function x=delinearize_around(obj,d_x,m_x)
            x = m_x  + d_x; 
        end
        
        % build a model from amount of resource types to utilization 
        function res=build_model(obj, inst_count_vec, util_vec)           
            n_st=1;% 1 state variable
            n_in=11;% 8 inputs
            n_out=1;% 1 output
            n_env=1; %environmental disturbances
            A=ones(n_st,n_st);
            B=ones(n_st,n_in);
            K=ones(n_st,n_env);
            C=ones(n_out,n_st);
            D=ones(n_out,n_in);
            
            [m_util_vec,              d_util_vec]= obj.linearize(util_vec);
            [m_inst_count_vec,  d_inst_count_vec] = obj.linearize(inst_count_vec); 

            u =  d_inst_count_vec;
            y = d_util_vec;

            data=iddata(y, u , obj.samplingTime); 
            m=idss(A,B,C,D,K,'Ts', obj.samplingTime,'DisturbanceModel','Estimate','InitialState','Estimate');
            m1=init(m);
            mm=pem(data,m1);
            
            res=struct( 'm',mm,... 
                              'm_inst_count_vec', m_inst_count_vec,...
                              'm_util_vec',m_util_vec); 
        end
        
        function yp=runModel(obj, res, inst_count_vec)
            % linearize input based on model steady state (linearxation
            % point) 
             d_inst_count_vec = obj.linearize_around(inst_count_vec, res.m_inst_count_vec); 
             data=iddata([],  d_inst_count_vec ,obj.samplingTime); 
             d_yp = sim(res.m,data,'noise');
             yp = obj.delinearize_around( d_yp.OutputData , res.m_util_vec);
        end
        
          function yp=predictWithModel(obj, res, inst_count_vec, util_interval )
              inst_count_interval = inst_count_vec(size(inst_count_vec,1)-120:end,:);
              d_inst_count_interval = obj.linearize_around(inst_count_interval, res.m_inst_count_vec); 
              d_util_interval = obj.linearize_around(util_interval, res.m_util_vec); 
             
             
             data=iddata(d_util_interval(1:60,:),  d_inst_count_interval(1:60,:) ,obj.samplingTime); 
             [e1,x0] = pe(res.m,data,'m(odel)'); 
             [e2,x0] = pe(res.m,data,'e(stimate)'); 
%              disp(strcat('improved error cov from ', num2str(cov(e1.OutputData)),...
%                  ' to ', num2str(cov(e2.OutputData))));
             res.m.x0=x0;
             data=iddata([],  d_inst_count_interval ,obj.samplingTime); 
             d_yp = sim(res.m,data,'noise');
             
             % d_yp = predict(res.m, data, 60)
             
             yp = obj.delinearize_around( d_yp.OutputData , res.m_util_vec);
          end
        
        % yp = predict(m,data,60)
        
        function instance_quantity_to_fire=get_policy_output(obj,policy,t)
            instance_quantity_to_fire = policy(t);
        end
        
        function res=plan(obj,system,duration,policy)          
                qty=[]; utils=[]; 
                for t=1:duration-1               
                    % disp(strcat('t=',num2str(t))); 
                    instance_quantity_to_fire = obj.get_policy_output(policy,t);
                    qty = [qty; instance_quantity_to_fire];
                    %utils = [utils;  system.getNextUtilAcc(instance_quantity_to_fire)];      %one server                    
                     system.getNextUtilAcc(instance_quantity_to_fire);
                end            
                res=struct('qty',qty);
        end
        
        function plan=test_plan(obj)
           
           plan=...
             [1, 1, 220
             60, 2, 120
             110, 3, 30
             160, 4, 10
             210, 5, 10
             260, 6, 10
             310, 7, 4
             360, 8, 15
             410, 9, 4
             460, 10, 3
             510,11, 5];
        end
        
        function plan=probe_plan(obj)
            plan=[];
%             for(i=)
%             ceil(rand*11)
        end
        
        function instance_quantity_to_fire=testPolicy(obj,t)
            default=[0 0 0 0 0 0 0 0 0 0 0 ];             
            % times = [1,60,110, 160, 210];
            % min, type, num
            plan=obj.test_plan();
                     
%             inst(1,:)=[15  0 0 0 0 0 0 0 0 0 ];
%             inst(2,:)=[0 64 0 0 0 0 0 0 0 0 0 ];             
%             inst(3,:)=[0 0 0  0 0 0 0 0 0 0 ];             
            
           ind=find(plan(:,1)==t);
           if isempty(ind)
               instance_quantity_to_fire=default;
           else
               % instance_quantity_to_fire=inst(ind,:);
               instance_quantity_to_fire = zeros(1,11);
               instance_quantity_to_fire(plan(ind,2))=plan(ind,3);
           end
        end
             
        function plot_instances(obj,system) 
            u=UtilityLib();              
            color_num =1;   
            instances = system.past_inst();
            % instances=instances';
            for k= 1 : size(instances,2)            
                plot(instances(:,k), strcat(u.lnstyle{mod(color_num,4)+1},u.color{mod(color_num,7)+1}),...
                                'LineWidth',2,...                             
                                'MarkerSize',10);       %  'MarkerEdgeColor','k',...          %    'MarkerFaceColor','g',...
                hold on;
                color_num = color_num+1;
            end
        end
        
        function acc_inst=accumulate(obj,acc_inst,t,instance_quantity_to_fire)
                    a=instance_quantity_to_fire;
                    if (t+60>size(acc_inst,2))
                       acc_inst=[acc_inst ; zeros(t+60-size(acc_inst,1),size(instance_quantity_to_fire,2))];          
                    end                   
                   acc_inst(t+1:t+60,:)=acc_inst(t+1:t+60,:)+ones(60,1)*instance_quantity_to_fire;
        end

        function res=predict_choice_outcome(obj,m, choice)
                system=obj.system_;
                % predict with moel
                current_time = system.t-1; 
                inst= system.past_future_inst();     
                inst_interval = obj.accumulate(inst( ( current_time-60:size(inst,1)) ,:),...
                                                          60,... %  current_time, 
                                                          choice); 

                past_util = system.past_util();
                util_interval =  past_util(((current_time-60):end),:);                 
                util_interval = [util_interval  ;  zeros(60,1)];                

                %  yp=obj.runModel( m, inst_interval);
                 yp=obj.predictWithModel( m, inst_interval,util_interval);                 
                 res_util = yp(size(yp,1)-60:end);
                 res_util_est = yp(size(yp,1)-120:size(yp,1)-60);                                   
                 
                 res=struct('inst',inst_interval,'util',res_util, 'util_est',res_util_est);
        end
                 
        function choices=generate_choices(obj)
            % a-star
            % ceil(240 ./ (prices_amz().get_compute_units))
            possibilities = ceil(240 ./ (prices_amz().get_compute_units));
            % 240   240    60    30    37    19    10    48    12     8      8
            A = {(0:1:40) (0:1:20) (0:2:60) (0:2:30) (0:2:37) (0:2:19) (0:1:10) (0:3:48) (0:2:12) (0:1:8) (0:1:8)};
            choices=    [ A{1}' zeros(length(A{1}),11-1)]; 

%             choices=vertcat(... %            [ A{1}' zeros(length(A{1}),11-1)],...
%             [ zeros(length(A{2}),1) A{2}' zeros(length(A{2}),11-2)],...
%             [ zeros(length(A{3}),2) A{3}' zeros(length(A{3}),11-3)],...
%             [ zeros(length(A{4}),3) A{4}' zeros(length(A{4}),11-4)],...
%             [ zeros(length(A{5}),4) A{5}' zeros(length(A{5}),11-5)],...
%             [ zeros(length(A{6}),5) A{6}' zeros(length(A{6}),11-6)],...
%             [ zeros(length(A{7}),6) A{7}' zeros(length(A{7}),11-7)],...
%             [ zeros(length(A{8}),7) A{8}' zeros(length(A{8}),11-8)],...
%             [ zeros(length(A{9}),8) A{9}' zeros(length(A{9}),11-9)],...
%             [ zeros(length(A{10}),9) A{10}' zeros(length(A{10}),11-10)],...
%             [ zeros(length(A{11}),10) A{11}' ]);
            
             %     [ (0:5:240)' zeros(size(0:5:240,2),10)],...            
            % choices = [ A{1}' zeros(length(A{1}),11-1)];
        end          
        
        function inst_cost = inst_cost(obj,inst_vec)
            inst_cost = inst_vec * (prices_amz().get_hourly_price) ;
        end
            
        function [cost,cost1,penalty]=cost(obj, inst_vec, util, inst_interval)              
              target_util = 60; % 0.8 utilizattion 
              util_violation =  util - target_util .* ones(size(util,1),size(util,2));
%               tmp(tmp>0)=0;
%               penalty = sum(-1*tmp) ; 
              util_violation(util_violation<0)=0;
              
              % hack: when model does not recognize there is no instance to
              % cover within the next control interval 
              for tt=1:obj.control_interval
                  if inst_interval(tt+60,:)==zeros(1,11)                      
                        util_violation(tt)=+inf;
                  end
              end
                  
              penalty = sum(util_violation) ; 
              cost1 =  obj.inst_cost(inst_vec); 
              % here we are adding (delta of dollars(instance hourly price
              % (e.g. 0.08 for small))) to (utilization deviation (e.g. 65))
              cost=100*cost1+1*penalty;              
%               disp(strcat('util:',num2str(util'))); 
%               disp(strcat('tmp:',num2str(tmp'))); 
%               disp(strcat('penalty:',num2str(penalty))); 
%               disp(strcat('cost1:',num2str(cost1))); 
%               disp(strcat('cost:',num2str(cost))); 

        end
        
        function instance_quantity_to_fire=mpcPolicy(obj,m,t)
            default=[0 0 0 0 0 0 0 0 0 0 0 ];             
            if (mod(t, obj.control_interval) ~= 0)
               instance_quantity_to_fire = default;
               return;
           end
           
           choices = obj.generate_choices();
           choice_opt =[];
           cost_opt=+inf; cost_inst_opt=+inf; penalty_opt=+inf; 
           for i=1:size(choices,1)               
                      choice=choices(i,:); 
                      res = obj.predict_choice_outcome(m,choice);       
                      [cost,cost_inst,penalty]=obj.cost(choice,res.util,res.inst);   % res.util is 60 values [system.t:system.t+60] 
                      if cost<cost_opt 
                          choice_opt=choice; 
                          cost_opt=cost;
                          cost_inst_opt=cost_inst;
                          penalty_opt=penalty; 
                      end
           end
           instance_quantity_to_fire = choice_opt; 
           disp(strcat('choice_opt:' , num2str(choice_opt))); 
           disp(strcat('cost_opt:' , num2str(cost_opt))); 
           disp(strcat('cost_inst_opt:' , num2str(cost_inst_opt))); 
           disp(strcat('penalty_opt:' , num2str(penalty_opt))); 
        end
        
       function test2(obj)
                obj.system_=model();
                system=obj.system_;
                experiment_duration=23*60
                training_duration = 9*60; 
                workload = system.initialize(experiment_duration); 
                res = obj.plan(system, training_duration,@obj.testPolicy);                  
                m=obj.build_model( system.past_inst(), system.past_util())                

                res=obj.plan(system,experiment_duration-training_duration, @(t) obj.mpcPolicy(m,t));    % sqr = @(x) x.^2;              
                
                tot_fired_inst  = sum(system.fired_inst,1);
                disp(strcat('tot_fired_inst:',num2str(tot_fired_inst)));
                tot_cost=obj.inst_cost(tot_fired_inst );
                disp(strcat('tot_cost:',num2str(tot_cost)));                
                disp('system.fired_inst are ...'); 
                for i=1:size(system.fired_t,1)
                   disp(strcat(num2str(system.fired_t(i)),'-->', num2str(system.fired_inst(i,:))));
                end
                
                subplot(2,1,1); plot(system.past_util()); 
                subplot(4,1,3); plot( system.past_future_inst()  ); 
       end
       
       function test1(obj)
                obj.system_=model();
                system=obj.system_;
                experiment_duration=540+60;                
                workload = system.initialize(experiment_duration); 
                res = obj.plan(system, 540,@obj.testPolicy);                  
                m=obj.build_model( system.past_inst(), system.past_util())
                
               past_util = system.past_util();  
               util_interval =  past_util(((system.t-1-60):end),:);                 
               
               u=UtilityLib();
               figure;
               subplot(2,1,1);   plot((system.t-60:system.t)', util_interval, u.color{1}); hold on; 
               inst =  system.past_future_inst();                   
               subplot(2,1,2); plot((system.t-60:system.t+60-1)', inst(system.t-60:system.t+60-1,:));
                                                                        
               choice=[170 0 0 0 0 0 0 0 0 0 0 ]; 
               res = obj.predict_choice_outcome(m,choice);
               subplot(2,1,1); plot((system.t:system.t+60)', res.util,  u.color{2}); 
               % subplot(2,1,2); plot((system.t-60:system.t+60)',res.inst);                
                % k = waitforbuttonpress 

              choice=[0 0 24 0 0 0 0 0 0 0 0 ]; 
               res = obj.predict_choice_outcome(m,choice);
               subplot(2,1,1); plot((system.t:system.t+60)', res.util,  u.color{3}); 
              % subplot(2,1,2); plot((system.t-60:system.t+60)',res.inst);                
               % k = waitforbuttonpress 

              choice=[0 0 0 0 0 0 0 0 0 1 0 ]; 
               res = obj.predict_choice_outcome(m,choice);
               subplot(2,1,1); plot((system.t:system.t+60)', res.util,  u.color{4});
              % subplot(2,1,2); plot((system.t-60:system.t+60)',res.inst);                

               subplot(2,1,1); line([system.t system.t],[0 100],...
                        'Color',u.color{1},...
                        'LineStyle',u.lnstyle{1},...
                        'LineWidth',2,...
                        'Marker','s',...
                        'MarkerEdgeColor','k',...
                        'MarkerFaceColor','g',...
                        'MarkerSize',6);                     
               
               % plot((system.t-60:system.t)', res.util_est);
                
%               subplot(2,1,2); plot((system.t-60:system.t+60)',res.inst);                
%                 subplot(4,1,3); plot( system.past_future_inst()  ); 
%                 subplot(4,1,4); plot(past_util); hold on; plot((system.t-60:system.t+60)',yp);  
       end      
        
       function test(obj)
                system=model();
                experiment_duration=540+60;
                
                workload = system.initialize(experiment_duration); 

                subplot(4,1,1);   plot(workload);           
                % apply test policy for duration=540
                res = obj.plan(system, 540,@obj.testPolicy);                  
                 subplot(4,1,2);  plot(system.past_util(),'g--');          
                 subplot(4,1,3);  plot(system.past_pu(),'g--');          
                 subplot(4,1,4);    obj.plot_instances(system);
                % up_to_now_acc_inst = system.past_inst();    %system.acc_inst(:,1:size(res.utils,1))';
                m=obj.build_model( system.past_inst(), system.past_util())
%                 
%                 choice=[60 0 0 0 0 0 0 0 0 0 0 ]; 
%                 u_future = obj.accumulate(system.past_future_inst(), system.t-1, choice); 
%                % u_future = system.past_future_inst(); 
%                 yp=obj.runModel( m, u_future)
%                     
%                 subplot(3,1,1);              
%                 past_util =  system.past_util();
%                 system.t
%                 plot((system.t-60:system.t-1)', past_util(system.t-60:end)); hold on;
%                 plot((system.t:system.t+60-1)', yp(system.t:system.t+60-1));
%                 
% %                subplot(3,1,2); plot(u_future)
% %                subplot(3,1,3); plot(yp);

        end      
        
%        The difference between prediction and simulation is that in prediction, the past values of outputs used for calculation are measured values while in simulation the outputs are themselves a result of calculation using inputs and initial conditions.
% The process of computing simulated and predicted responses over a time range starts by using the initial conditions to compute the first few output values.                 
%         model = update_identified_model
%         goal_path=indentify_trajectory(current_val)??
%         current_action=argmin(|goal_path-model|+b|cost|)
%         perform_action
%         
    end
    
end

